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COMPUTER SIMULATION OF A CYLINDRICAL GALAXY 


By Frank Hohl 
Langley Research Center 

SUMMARY 

A two-dimensional computer model is used to investigate the evolution of a self- 
gravitating cylindrical galaxy. The evolution of systems containing up to 100 000 mass 
rods is investigated by following the self-consistent motion of the mass rods. As the 
system evolves in time, it displays the various shapes observed in actual galaxies; how- 
ever, any spiral structure which appears during the initial stages of the evolution disap- 
pears after about 10 rotations, and the final state of the system is that of a rotating bar 
or a circular cluster with a high central condensation. It is found that violent relaxation 
or phase mixing causes any unstable system quickly to build up large random velocities. 
The rate of buildup of random motions in stable systems was found to be described by a 
Chandrasekhar-type relaxation time. 

INTRODUCTION 

Two-dimensional computer models have recently been used by Hockney (ref. 1) and 
by Hohl (refs. 2 to 4) to simulate the evolution of collisionless self-gravitating systems. 
In the two-dimensional model, the N stars in the system are represented by infinitely 
long mass rods which have a mass m per unit length. These stars are then allowed to 
move in the two directions perpendicular to their axes. This approximation leads to a 
completely two-dimensional problem and greatly simplifies the calculation of the gravi- 
tational potential. 

In the present calculations the number of stars in the system has been greatly 
increased as compared to the calculations previously made (ref. 4). The results indi- 
cated that especially for systems containing small numbers of stars, the detailed struc- 
ture displayed by a system depended markedly on the number of stars in the system. 

Also, the systems quickly build up temperature or random motion of the stars which 
destroys any spiral structure. Tests were performed to determine whether the randomi- 
zation was caused by collisional effects introduced by the model or by the violent relaxa- 
tion proposed by Lynden-Bell (ref. 5) and investigated by Hohl and Campbell (ref. 6). 


SYMBOLS 


major axis; also, acceleration 

minor axis 

distance 

gravitational field 
gravitational constant 
mass per unit length 
total mass per unit length 
star density 
total number of stars 
potential energy 


time 

kinetic energy 
total energy, T + P 
velocity 

rms value of relative velocity in balanced cylinder 
transverse velocity 

positions along x and y axes, respectively 


radius, 



complex coordinate, x + iy 


y convergence parameter 

p mass density 

r c relaxation time 

Tg equilibrium rotation period 

cp gravitational potential 

oj frequency of rotation 

o>g equilibrium frequency of rotation 

Subscripts: 

i,j,k,n,m summation indices 

max maximum 

min minimum 

r radial 

x,y x- and y- component 

6 azimuthal 

/ ) denotes expectation values 

MODEL AND METHOD OF SOLUTION 

The model consists of a large number of mass rods. The position and velocity of 
all the mass rods is stored in the computer. Newton's laws of motion and Newton's 
theory of gravitation are used to advance the position and velocity of all the stars step- 
wise in time. 

The area which contains the galaxy of stars is divided into a square array of cells 
(usually 101 x 101). At the center of each cell a mass density p n m is defined which is 
given by the number of stars in that cell. The integers n and m increase in the 
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x- and y-direction, respectively. The gravitational potential <p n>m is then obtained 
from the density p n m by solving the two-dimensional Poisson equation 

= 47rGp(x,y) (1) 

8x^ 9y2 

by finite difference methods (ref. 4). The standard five-point difference equation (ref. 7) 
^n+l.m + ^n,m+l + ^n-l,m + ^m-l " 4<p n,m = 47rGp n,m ^ 


was used to solve for the potential distribution. The cell dimensions Ax and Ay are 
taken to be equal to unity. This set of simultaneous equations is solved by an iteration 
method on the Control Data 6600 computer in the form 


„r+l 


^n,m ^n,m + 




m 


+ 


n+l,m + ^n,m-l " r ^n,m+l 


r+1 


+ <Pr 


- 4 cp - 


n,m 


4vGp n> 


m 


( 3 ) 


For the purpose of saving computer storage and increasing the convergence rate, the new 
values of cp (that is, <p r+ l) which have already been calculated during a particular 
iteration are used in the right-hand side of equation (3). The superscript r refers to 
the rth iteration and the parameter y is adjusted to give the maximum rate of conver- 
gence. The potential at the boundary of the rectangular region is required in the solution 
of the Poisson equation. At an arbitrary boundary point a distance z = x + iy from the 
center of the mesh, the potential is given by 


<p(z) - 2G 


p n,m l°^e 


z - z 


n,m 


= 2GM log e |z 


+ 2G 


n,m 


I 

n,m 


p n,m Re 


log. 


.(* - z -¥) 


( 4 ) 


where 

the cell 
as 


M is the total mass in the system and z n>m = x 

z ixi I 

n,m. Since P n>m is nonzero only for — ^ — 


n,m 

< 1 , 


+ iy„ is the coordinate of 
J n,m 

equation (4) can be written 


where 


cp(z) = 2GM log e |z| - 2GRe 


I 


a k 

kzk 


( 5 ) 


and the series expansion for 


l°g e 


a k ^ p n,m z n,m 
n,m 

^1 — l* as l> een truncated after 15 terms. 


( 6 ) 
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If the motion of all the stars in the system is advanced for a small time step 6t, 
the mass distribution p n m will not change very much. The change in the gravitational 
potential will then also be very small. Thus, the solution of the finite difference form of 
the Poisson equation (eq. (3)) by an iteration method which uses the potential from the pre- 
vious cycle as an initial guess will converge very rapidly. The accuracy of the iterative 
solution of the Poisson equation is continually checked during the calculations. This veri- 
fication is made by obtaining the values of the potential and the field at several selected 
points by summing directly the contribution from each star. The values so obtained 
agree with those obtained from the solution of the Poisson equation to at least the first 
three digits. The number of iterations required for a 51 x 51 mesh was 5 to 7 and for a 
101 X 101 mesh, 12 iterations. The calculation time for one cycle of a 10 000-star sys- 
tem with a 101 X 101 mesh is 7 seconds on the Control Data 6600 computer. This time 
includes such operations as checking of the gravitational potential, calculation of the 
energy and angular momentum, and writing the positions and velocities of all stars on 
tape. 

The initial guess of the potential at t = 0 was determined by using the analytical 
expression for the potential of a homogeneous elliptical cylinder. (See ref. 8.) If the 
cylinder containing N stars of unit mass is described by 

x 2 y2 

^r+ i-= 1 
a^ b^ 


the potential interior to this cylinder is given by 


< p (*, y) = + tt) + <p( 0 >°) 


a + b \a b 


where 


<p(0,0) = 2G ^ P n ,m lo e e | z n, 


m 


n.m 


and z n m is the position of a cell measured from the center of the mesh. For points 
exterior to the cylinder, the potential is 


<p(x,y) = 2GN 


a ^ + k + \j{ b^ + k)( a^ + k) b^ + k + \j(b^ + k)( a^ + k) 


where 




^2 + y2 _ a 2 - b^ + ^(a2 + b^ - x^ - y^j - 4^a^b^ - b^x^ - a^y2) 
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The motion of the stars is described by the differential equations 


dV ? 

dt 

dV 

~dt" 


9 cp 
dx 


y _ 


and 


9 cp 

’ 8 y 


cbT 




( 7 ) 


Vx = 5F 


v 

v y _ dt^ 




( 8 ) 


For a star in the (n,m)th cell equations (7) and (8) are in finite difference form 

V x( t + f)' 


_1 

6t 


and 


2( < ’ C n-l,m ” ^n+l,m) 

(9) 

■ * 1 ) 

(10) 


with similar equations for the y-components. As given by equation (9), all stars in the 
cell n,m experience the same gravitational field. For the calculations reported in the 
present paper, the field acting on a star was obtained by means of a bilinear interpolation 
of the fields at the four cell centers surrounding the position of the particle. The right- 
hand side of equation (9) is then replaced by the interpolated field. The new velocity 
V (t + and position x(t + 6t) are then obtained from those at the beginning of the 
time step by means of the equations 

Vx(t + f) = V x (t - f) + f [(1 + 6y)(l - 6x)(«» n+ljm - <?„_!,„) 

+ 5y(l - 6x)(v n+l m+1 - Vl,m+l) + 6x(1 ‘ 6y) CW,m ‘ Vm) 

+ 6 x 5 y (® n+2 , m+1 + y’n.m+l)] 


»1) 


and 


x(t + 6t) = x(t) + 5tV x (t + y) 


( 12 ) 
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One complete cycle for advancing the motion of the system by a time 6t consists 
of the following steps. First, the distribution of mass p n m is used to obtain the gravi- 
tational potential cp n m by numerically solving the Poisson equation. Second, the grav- 
itational field at the position of the stars is computed from the potential <p n)in . Third, 
Newton's laws of motion are used to advance the motion of all the mass rods for a small 
time step 6t. These three steps represent one cycle and are repeated until the desired 
evolution of the system is achieved. 

The kinetic energy of the system at time t is 



where the summation extends over all stars and the potential energy is 

P = II W (14) 

n,m 

The total energy U of the system is given by 

U = T + P (15) 

For most of the systems investigated the calculated total energy during a given 
computer run was conserved to within 0.5 percent and the angular momentum remained 
practically constant. 


RESULTS AND DISCUSSION 


Computer simulations have been performed for a number of systems with an ini- 
tially uniform distribution over a circular or elliptical region in x,y- coordinate space, 
zero thermal velocity, and various values of initial solid-body rotation. In particular, 
systems with a ratio of major axis a to minor axis b of 3 and with an initial solid- 
body rotation given by 


o o 



8 MG 
(a + b) 2 


(16) 


have been investigated in detail. The o>g given by equation (16) will initially approxi- 
mately balance the gravitational attraction. It was found that for a system containing 
2000 stars with a 51 x 51 grid to solve the Poisson equation, randomization effects 
destroyed the spiral structure after only 2 to 3 rotations. The results are shown in fig- 
ures 1 to 3. The normalization G = m = 1 was used in the calculations and the time 
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is shown in rotation periods 


r =^H- 
g “ g • 


Figure 1 shows the evolution of the 2000- star 


system in x,y-coordinate space. These results show that any structure which appears 
initially quickly randomizes after less than 3 rotations. The evolution of the velocity 
distribution in V x ,Vy- space is shown in figure 2. Figure 3 shows the evolution of the 
velocity distribution in coordinates that rotate with the equilibrium angular velocity a> g . 
The velocity (V q - rco g ) where Vq is the azimuthal velocity and r is the radius from 
the center of the system to a star is plotted as a function of the radial velocity V r . The 
velocities in figures 2 and 3 are plotted on nearly equal scales. Initially, the thermal 
(or random) velocity of all mass rods is zero since initially the system has a solid-body 
rotation only. As the system evolves in time the thermal velocity builds up quickly and 
the system expands in (V q - ro>g),V r -space. After a time t = 1.75Tg the random veloc- 
ities are of the same order as the initial circular velocities and there is little further 
change in the velocity distribution. The long tail in the velocity distribution in figure 3 
is due to differential rotation. These results indicate that for the calculations reported 
by Hockney (ref. 1) with 2000 stars and a 48 x 48 grid, the system randomizes in about 
2 rotations. The same holds for the results given by Hohl (refs. 3 and 4) for 2000 stars 
(with a 51 x 51 grid). To reduce the randomization effects (i.e., fluctuations of order 
l/^N)the calculations were repeated with an identical initial elliptical distribution but with 
4000 stars and a 101 x 101 grid to solve the Poisson equation. The results are shown in 
figures 4 and 5. Figure 4 shows that the spiral structure now has a somewhat greater 
stability but still tends to randomize after about 4 rotations. The corresponding evolu- 
tion in V x ,Vy-velocity space is shown in figure 5. The calculations were repeated with 
the number of stars increased to 10 000; the resulting evolution in x,y-coordinate space 
is shown in figure 6. Figure 6(b) shows that at t = 6r g the spiral structure begins to 
randomize. From these results it appears that in order to simulate collisionless sys- 
tems for the first n rotations, nearly 2n x 10^ stars are required. 


An interesting point is that all systems investigated which have an initially elon- 
gated distribution and a given solid-body rotation go through an initial evolution similar 
to that shown in figure 6(a) up to t = 1.0T g . At t = 0.75T g the system shows a struc- 
ture nearly identical to that of the peculiar galaxy (GB 1) discussed by Burbidge, Burbidge, 
and Shelton (ref. 9). Another such integral- sign- shaped galaxy is NGC 4656/7 reproduced 
in the Hubble Atlas (ref. 10, p. 40). Figure 6 shows that at t = 3.0r g the two spiral 
arms begin to overlap and give the appearance of a circular ring of stars separate from 
the nucleus, presenting a shape similar to NGC 3145 (ref. 10, p. 21). At t = 4.0jg the 
system shows a structure similar to NGC 5194 (ref. 10, p. 31). Furthermore, at 
t = 3.0Tg the spiral structure has almost disappeared and at t = 3.75Tg, after less than 
one rotation, the spiral structure is again very pronounced. This phenomenon appears 
to reoccur from t = 4.50r g to t = 6.00T g . However, at t = 6.00r g randomization 
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effects have already appreciably affected the evolution of the system. This behavior 
indicates that the spiral arms may be short lived but re-form repeatedly as in the theory 
of Goldreich and Lynden-Bell (ref. 11). Since the spiral arms do not remain for a long 
time it is difficult to determine whether the stars that are in a particular arm remain 
there, as in the theory of Goldreich and Lynden-Bell, or whether they move from arm to 
arm as would be expected from the theory of Lin and Shu (refs. 12 and 13). The evolu- 
tion of the velocity distribution for the 10 000-star system is shown in figure 7. At 
t = 0.25rg the system shows an interesting structure. Condensation in velocity space 
has occurred and striation appears in the velocity distribution. The corresponding evo- 
lution in (Vg - ru>g),V r -velocity space is shown in figure 8. This distribution shows little 
change after the first 3 rotations. Again, the scales for plotting the velocity distributions 
in figures 7 and 8 are about the same. 

The number of stars in the system was finally increased to 100 000 to determine 
whether a more stable spiral structure would result. The evolution in x,y- coordinate 
space for the 100 000-star system is shown in figure 9. These results show that no 
stable spiral structure develops even for large numbers of stars in the system. 

The effect of giving the system an initial rotation 30 percent larger than that given 
by equation (16) is shown in figure 10. There are 10 000 stars in the system. Figure 10 
shows that the initially increased speed of rotation causes the system to break up into 
clusters and evolve more irregularly as compared with the system shown in figure 6. 
After 15 rotations the system has a bar shape with a diffuse cloud of stars around the 
bar. The long-time evolution of various systems investigated shows that the bar shape 
structure is favored as the final state of the system. 

The effect of placing a large mass at the center of the system was also investigated. 
Figure 11 shows the positions of the 10 000 stars in the system at different times. The 
central mass has 10 times the mass of the 10 000 stars in the system. The initial dis- 
tribution was chosen to produce a spiral structure. The time is given in arbitrary units. 
The figure shows that the spiral pattern simply winds up because of differential rotation. 
No persistent spiral pattern is found. 

The evolution and the final state of a system might be expected to depend on the 
initially elongated distribution. The evolution of systems with initially circular distribu- 
tions in coordinate space were therefore investigated. Figure 12 shows the evolution of 
a 10 000- star system with an initially uniform distribution over a circular region and 
with an initial solid-body rotation equal to that required to balance gravitation. The 
gravitational field g inside such a system is 

g = 27rGpr (17) 

To balance this gravitational field requires a centrifugal force equal to u^r where 
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a) = \l2nGp 


(18) 


The rotational period for such a system is therefore 

_ 2tt 

Tg “ \Z2 ttGp 


(19) 


Figure 12 shows that for the first 2 rotations the evolution of the balanced cylinder 
containing 10 000 stars is very similar to that of the 2000-star system investigated by 
Hockney (ref. 1) and by Hohl and Park (ref. 4). For these 2000- star systems of refer- 
ences 1 and 4, the results showed that the relatively small fourth-mode surface wave 
around the periphery of the cylinder quickly changed to a second- mode surface wave 
(egg shape), then all structure was lost, and the system finally assumed a circular shape. 
This sequence is due to the large randomizations (pressure) that developed to balance the 
gravitational attraction towards the center of the system. Figure 12 shows that the evo- 
lution is quite different when randomization effects are reduced by increasing the number 
of particles to 10 000. The fourth-mode surface wave now continues to grow and after 
5 rotations a four- arm spiral begins to form. After about 8 rotations two of the arms 
have disappeared and the system has assumed a bar shape. Figures 13 and 14 are plotted 
on the same scale and they show that after about 9 rotations the random or thermal veloc- 
ities have increased to such an extent that they are of the same magnitude as the rota- 
tional velocities. The evolution of random velocities shown in figure 14 shows an 
almost linear increase in the random or thermal motion of the mass rods. Figure 12 
shows that the second-mode surface wave will finally dominate and the system assumes 
a bar shape which appears to be the most probable state for the system. Figure 15 shows, 
again, the evolution of a 10 000 -star system like that shown in figure 12. The only dif- 
ference is that a different pseudo-random number sequence was used for the initial posi- 
tions of the stars. The evolution is very similar to that shown in figure 12 and again the 
final state is a bar-shape system. When the trajectories of individual stars are plotted 
they are found to oscillate in the potential well set up by the high mass concentration in 
the bar. 


It is of interest to determine the extent to which randomization or collisional effects 
are present in the computer model. The relaxation time tc (ref. 14) can be approxi- 
mated by applying the method for point stars (ref. 15) to mass rods. 


Consider an encounter between a rod of mass m 2 and velocity V with a station- 
ary rod of mass m^. If D is the distance of closest approach (impact parameter) then 
the transverse force acting on the moving star m 2 during a time 2D/V can be approx- 
imated as 2Gm 1 m2 / /D. The moving star has therefore acquired a transverse velocity 


6V-p = 


4Gmj 

V 


( 20 ) 
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For a star interacting with a system containing many stars, the effect of many individual 
encounters must be summed. The number of encounters in a time t with impact param- 
eter between D and D + dD is tVn dD where n is the density of stars in the sys- 
tem. The expectation value of V,p 2 is then given by 


< V T 2 ) = 


'max 


'min 


16tnG 2 mi 2 _ 16tnG2 mi 2(D ma x “ D min )_ lGtnG^D^ 
V V V 


( 21 ) 


since in general D max » D m ^ n where D max is the dimension of the system and 
D m in I s the cell size used in the calculations. The relaxation time r c is defined to 
be the time required for to be of the same order of magnitude as V 2 . Equa- 

tion (21) then becomes 

Vg 3 

t c = 1 (22) 

16pmG D max 

where Vg is the root- mean- square value of the velocity of the stars and D max is the 
radius of the system. Consider now a uniformly rotating cylinder of stars of radius 
D m ax and constant density p. The period of rotation is given by equation (19) and the 
mean quadratic difference of velocity between two stars is 

Vg = D max (27rGp) t (23) 


The ratio of the relaxation time r c to the rotational period Tg is then 

2 

T c _ ff ^max n _ N / 04 \ 

T g 8 " 8 ( 

Equation (24) is a rough estimate but it indicates that for a system containing 
10 000 stars, randomization effects would be expected to occur after about 1000 rota- 
tions. Since the observed randomization for the unstable system occurs after only a few 
rotations, the Chandrasekhar-type (ref. 14) calculation (eq. (22)) of the relaxation time 
does not apply. For a system with an initial solid-body rotation there is only one veloc- 
ity at any given point. Thus, in the neighborhood of any particular star very little rela- 
tive motion is present and the usual assumption of random motions underlying the calcu- 
lation of the relaxation time does not apply. It was suggested by Dr. M. Henon of 
Observatoire de Nice (France) in a personal communication to the author that the fol- 
lowing estimate of the relaxation time may be more appropriate. The acceleration at a 
distance D away from a single star is 2Gm/D. At the center of the system the expec- 
tation value of the square of the total acceleration a is therefore 
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fD 


<a2> = 


max 


4G 2 


D 


min 


D 2 


m2 2-nriD dD = 8irpG 2 m log e ° max 


D 


mm 


(25) 


If the acceleration is taken to be constant, the perturbations in the velocity of the stars 
<V t 2 ^> will grow as the square of the time 

<V T 2 ) = t 2 (a2> 

This result is in agreement with the results shown in figure 14. In the usual case 
grows linearly with time. The relaxation time r c is defined to be the time required 
for <V'p 2 y’ to be of the same order of magnitude as Vg 2 . Thus, 



and 


N 


,1/2 


T g 4w l 


log 


D 
e D 


max 


(26) 


miny 


The distance D m j n is equal to the dimension of a cell so that max 

l^min 

T c_ # 
r g = 25 


51 and 

(27) 


For N = 2000, 4000, and 10 000 the corresponding values of r c /Tg are 1.8, 2.5, and 
4.0, which are in good agreement with the numerical results. This type of fast relaxa- 
tion which occurs in unstable systems is the violent relaxation or phase mixing described 
by Lynden-Bell (ref. 5). 


The results given in figures 12 and 15 show that the uniformly rotating cylinder is 
unstable. Figure 16 shows the effect on the stability when one-half of the stars rotate 
clockwise and the other half rotate counterclockwise. The cylinder is now stable and 
retains its original shape even after many rotations. Figure 17 shows the variation of 
the total kinetic energy for the stable cylinder shown in figure 16. The initial small 
oscillations are quickly damped out. The variation of the kinetic energy for the unstable 
systems shown in figures 12 and 15 is essentially identical to that shown in figure 17. 
Because the system of counterrotating stars is stable, it is suited for determining the 
rate at which thermalization effects occur in the two-dimensional model. These effects 
are calculated by determining the rate at which equipartition of energy occurs among 
stars of different masses. Initially, the balanced system contains a uniform distribution 
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Q 

of 5000 stars of mass m-t = — m rotating clockwise plus a uniform distribution of 

z 

5000 stars of mass m 2 = i m rotating counterclockwise. The resulting evolution of the 

z 

system is shown in figures 18 and 19. The normalized kinetic energy for the heavy stars 
is now defined as 

Tj* = ^ 1 7 Vi 2 (heavy stars) 

i 

where the summation extends over the 5000 heavy stars. Similarly for the light stars 

T 2* = 4 = iI V 3 2 (light stars) 

J i 

* * 

Figure 20 shows the percent variation of the difference between T^ and T 2 . Ini- 
tially, the heavy stars have an excess of kinetic energy because of their larger mass. 

The total unnormalized kinetic energy of the heavy stars is ^ and the unnormal- 

1 * * 
ized kinetic energy of the light stars is — mT 2 . The fast drop in the energy Tj 

during the first 2 rotations is probably due to collective effects, as are the oscillations 
shown after 4 rotations. However, figure 20 clearly shows that the heavy stars are 
losing kinetic energy, and after 11 rotations the average value of T 2 * is 0.5 percent of 
the total kinetic energy larger than Tj*. For equipartition of energy between the heavy 
and the light stars, T 2 * = 3Tj*. If the same rate of increase in T 2 * as shown in fig- 
ure 20 is assumed, then equipartition would be achieved after about 1100 rotations, which 
is in good agreement with the value 1250 given by equation (24) for N = 10 000. 

CONCLUDING REMARKS 

The results obtained with the two-dimensional model for a self-gravitating cylin- 
drical galaxy show that as the system evolves in time, it displays the various shapes 
observed in actual galaxies. However, the rod-star approximation is likely to be valid 
only for very few galaxies, such as the cigar-shaped galaxy NGC 2685. In order to simu- 
late the evolution of systems like spiral galaxies, calculations are now in progress for 
systems consisting of finite length mass rods and point masses moving in a plane. 

It was found that for stable systems the Chandrasekhar-type calculation of the 
relaxation time was applicable and that the stable two-dimensional model can indeed be 
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considered collisionless. The fast randomization which occurred for unstable systems 
must be ascribed to violent relaxation or phase mixing as predicted by Lynden-Bell. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., March 3, 1969, 

129-02-01-01-23. 
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Figure 1.- Evolution of a 2000- star system in coordinate space. 
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Figure 2.- Evolution of a 2000-star system in V x ,Vy-velocity space. 
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Figure 4.- Evolution of a 4000-star system in x,y-coordinate space. 
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Figure 4.- Concluded. 
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(a) Evolution up to t = 2.75ig. 

Figure 5.- Evolution of a 4000-star system in V x ,V y -ve locity space. 
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Figure 5.- Concluded. 
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(a) Evolution up to t = 2.75ig. 

Figure 6.- Evolution of a 10 000-star system in x,y-coordinate space. 
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(b) Evolution up to t = lO.OOig. 
Figure 6.- Concluded. 
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(a) Evolution up to t = 2.75tg. 

Figure 7.- Evolution of a 10 000-star system in V x ,Vy-velocity space. 
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(a) Evolution up to t = 2.75ig. 

Figure 8.- Evolution of a 10 000-star system in (V0 - velocity space. 
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(b) Evolution up to t = lO.OOig. 
Figure 8.- Concluded. 
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Figure 9.- Evolution of a 100 000-star system in x,y-coordinate space. 
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(b) Evolution up to t = 4.6tg. 
Figure 9.- Continued. 
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(c) Evolution up to t = 10.5ig. 
Figure 9.- Concluded. 
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(a) Evolution up to t = 1.37ig. 

Figure 10.- Evolution in x,y-coordinate space of a 10 000-star system with an initial rotation equal to 1.3wg. 
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(b) Evolution up to t = 5.5ig. 
Figure 10.- Continued. 
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Figure 10.- Concluded. 
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Figure 11.- Evolution of a 10 000-star system in the presence of a large central mass equal to 10 times the mass of the 10 000 stars. 
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Figure 11.- Concluded. 
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(a) Evolution up to t = 2.75ig. 

Figure 12.- Evolution of a 10 000-star system in x,y-coordinate space. The initial solid-body rotation of the system equaled that 

required to balance the gravitational attraction. 
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(b) Evolution up to t = 12.5ig. 
Figure 12.- Concluded. 
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(a) Evolution up to t - 2.75tg. 

Figure 13.- Evolution of a 10 000-star system in V x ,Vy-velocity space where initially the centrifugal force balances 

the gravitational attraction. 
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(b) Evolution up to t = 12.5ig. 
Figure 13.- Concluded. 


41 



(a) Evolution up to t = 2.75ig. 

Figure 14.- Evolution of a 10 000- star system in V r ,(V 0 - ro> g ) -velocity space. The initial solid body rotation of the system 
equaled that required to balance the gravitational attraction. 
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Figure 14.- Concluded. 
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Figure 15.- Evolution of a balanced system like that shown in figure 12, but with a different pseudo-random number sequence 

for the initial star positions. 
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Figure 16.- Evolution of a balanced system with 5000 stars rotating clockwise and 5000 stars rotating counterclockwise. 







45 


05 



Figure 17.- Evolution in time of the kinetic energy for the balanced cylinder shown in figure 16. 
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Figure 18.- Evolution in x,y-coordinate space for a stable system with counterrotating stars of different mass. 
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Figure 19.- Evolution of a stable system of counterrotating stars in V x ,Vy-velocity space. 
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Figure 20.- Rate at which equipartition of energy occurs for a stable system containing 5000 stars of mass m/2 and 

5000 stars of mass 3m/2. 
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